module ModuleDistributions

    use Globals
    use Functions_wa
    use Toolbox
    
contains
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution of tax-and-transfer rates, income, net worth
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine income_distribution

        implicit none

        ! Local variable declarations
        real(8), dimension(NS) :: yk, yl, ytot, tax_paid                                        ! capital, labor, and total income
        real(8), dimension(NS) :: ytot_sorted, yl_sorted, tax_paid_sorted, transfers_sorted     ! total income, labor income, tax paid, transfers (all sorted)
        real(8), dimension(NS) :: mu_sorted, mu_cumu, a_sorted, total_tax, total_tax_sorted     ! measure sorted, cum. measure, assets sorted, cap.+lab. tax, also sorted
        real(8), dimension(NS) :: av_tax_sorted, marg_tax_sorted, av_tax                        ! avg. tax sorted, marg. tax sorted, avg. tax 
        real(8), dimension(NS) :: marg_tax, dmarg_tax, dmarg_tax_sorted                         ! marg. tax, change in marg. tax, change in marg. tax sorted
        integer :: ind_sort(NS)                                                                 ! index for sorting
        integer :: Nstate                                                                       ! auxiliary variable equal to NS
        real(8), dimension(NS) :: muu, S1u, S2u                                                 ! auxiliary variables equal to measure, asset states, productivity states
        real(8) :: avg_tax_q(5), avg_transfers_q(5), inc_share(7), a_share(7)                   ! avg. tax-transfer by quintile, avg. transfer by quintile, inc. share by quintile and top 10 and top 1, asset shares
        real(8) :: avg_income_tax_q(5), a_share_top10, a_share_top1                             ! avg. lab. tax, asset shares top 10 and top 1
        real(8) :: avg_income_tax_ind_q(5), avg_transfer_ind_q(5)
        real(8) :: avg_cap_income_tax_q(5), avg_total_tax_q(5), avg_tax_ind_q(5)                ! avg. cap. tax, avg. total tax, avg. tax other concept
        real(8) :: marg_tax_ind_q(5), dmarg_tax_ind_q(5)                                        ! marg. tax by quintile, change in marg. tax by quintile
        real(8) :: inc_share_top10, inc_share_top1                                              ! inc. share top 10 and top 1
        integer :: p_20_ind, p_40_ind, p_60_ind, p_80_ind, p_90_ind, p_99_ind, p_50_ind         ! indices for percentiles
        integer :: p_1_ind, p_5_ind, p_10_ind                                                   ! indices for percentiles
        integer :: iter, is                                                                     ! counter
        real(8) :: chit                                                                         ! auxiliary variable transfer
        real(8), dimension(NS) :: et, transfers                                                 ! auxiliary variable transfer, transer
        real(8), dimension(NS) :: income_tax, income_tax_sorted                                 ! income tax, and sorted
        real(8), dimension(NS) :: capital_income_tax, capital_income_tax_sorted, yk_sorted      ! cap. tax, sorted, cap. income sorted
        real(8) :: auxx_pareto(NS,2), auxy_pareto(NS,1), coef_pareto(2,1), auxx_trans(2,NS), aux22(2,2), aux2NS(2,NS)
        real(8) :: yl_share(5), yk_share(5)                                                     ! labor and capital income shares by total income quintiles

        ! auxiliary variables
        Nstate = NS     ! number of individual states
        muu = mu        ! measure over individual states
        S1u = S(:,1)    ! assets for individual states
        S2u = S(:,2)    ! productivities for individual states

        ! Income variables
        yk = r*S1u          ! capital income
        yl = wge*S2u*h_pol  ! labor income
        ytot = yk+yl        ! total income
    
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Tax and transfer rates
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        ! Total taxes and transfers
        tax_paid = TAX(yl,yk)   ! total tax-and-transfer paid
        av_tax = tax_paid/ytot  ! total tax-and-transfer average rate

        do is = 1, Nstate
            marg_tax(is)  = dTAX_sc(yl(is),yk(is))      ! total tax-and-transfer marginal rate
            dmarg_tax(is) = ddTAX_sc(yl(is),yk(is))     ! total tax-and-transfer change in marginal rate
        enddo

        ! Transfers
        et   = exp(-rt*((ytot)/ymean-omegat))                   ! auxiliary variable
        chit = exp(rt*omegat)/(1.0D0 + exp(rt*omegat))          ! auxiliary variable
        transfers = (mt*ymean) * (et/(1.0D0+et)) * (1.0D0/chit) ! transfer
    
        ! Labor income tax
        income_tax = exp(log(lambda)*((yl/ymean)**(-gamma)))*yl ! labor income tax
    
        ! Capital tax
        capital_income_tax = tauk*yk                            ! capital income tax
    
        ! Total tax (without transfers)
        total_tax = income_tax + capital_income_tax             ! total tax (without transfers)

        ! Aggregate transfer payments
        Tagg = sum(transfers*muu)           ! aggregate transfers
        print *, ''
        print *, 'T/Y in percent =', Tagg/Yagg*100d0

        ! sort total income array
        ytot_sorted = ytot
        ind_sort = 0
        call hpsort_eps_epw(Nstate, ytot_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income
        tax_paid_sorted = tax_paid(ind_sort)                        ! total tax-and-transfer sorted by total income
        transfers_sorted = transfers(ind_sort)                      ! transfer sorted by total income
        income_tax_sorted = income_tax(ind_sort)                    ! labor tax sorted by total income
        capital_income_tax_sorted = capital_income_tax(ind_sort)    ! capital tax sorted by total income
        total_tax_sorted = total_tax(ind_sort)                      ! total tax (without transfers) sorted by total income
        av_tax_sorted = av_tax(ind_sort)                            ! average tax-and-transfer rate sorted by total income
        marg_tax_sorted = marg_tax(ind_sort)                        ! marginal tax-and-transfer rate sorted by total income
        dmarg_tax_sorted = dmarg_tax(ind_sort)                      ! change in marginal tax-and-transfer rate sorted by total income


        print *, 'minimum income in dollars = ', ytot_sorted(1)/ymean*90185D0, mu_sorted(1)
        print *, 'minimum income = ', ytot_sorted(1), wge, r, ymean
        print *, 'transfers received at the minimum = ', transfers_sorted(1)/ymean*90185D0
        print *, 'transfers net of taxes at the min (incl capital tax) = ', -tax_paid_sorted(1)/ymean*90185D0
        print *, 'transfers net of taxes at the min (excl capital tax) = ', (transfers_sorted(1)-income_tax_sorted(1))/ymean*90185D0
        print *, 'transfers in m terms  = ', mt

        print *, 'transfers received at the minimum = ', transfers_sorted(1)/ymean*90185D0*(ymean/0.791095232895847D0)
        print *, 'transfers net of taxes at the min (incl capital tax) = ', -tax_paid_sorted(1)/ymean*90185D0*(ymean/0.791095232895847D0)
        print *, 'transfers net of taxes at the min (excl capital tax) = ', (transfers_sorted(1)-income_tax_sorted(1))/ymean*90185D0*(ymean/0.791095232895847D0)
        print *, 'transfers at exactly zero income  = ', mt*ymean/ymean*90185D0*(ymean/0.791095232895847D0)
        print *, 'transfers in m terms  = ', mt
        print *, 'transfers in m terms, rescaled  = ', mt*(ymean/0.791095232895847D0)

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        ! "Quintiles"
        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        p_99_ind = minloc(abs(mu_cumu-0.99d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.50d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        ! compute average tax rates by quintile (CBO concept: total tax of the quintile divided by total income of the quintile)
        avg_tax_q(1) = sum(tax_paid_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_tax_q(2) = sum(tax_paid_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_tax_q(3) = sum(tax_paid_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_tax_q(4) = sum(tax_paid_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_tax_q(5) = sum(tax_paid_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        ! compute average tax rates, at the indiv level (mean of average tax-and-transfer rate)
        avg_tax_ind_q(1) = sum(av_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_tax_ind_q(2) = sum(av_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_tax_ind_q(3) = sum(av_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_tax_ind_q(4) = sum(av_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_tax_ind_q(5) = sum(av_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute marginal tax rates by quintile (mean of marginal tax-and-transfer rate)
        marg_tax_ind_q(1) = sum(marg_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        marg_tax_ind_q(2) = sum(marg_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        marg_tax_ind_q(3) = sum(marg_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        marg_tax_ind_q(4) = sum(marg_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        marg_tax_ind_q(5) = sum(marg_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute change in marginal tax rates by quintile
        dmarg_tax_ind_q(1) = sum(dmarg_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        dmarg_tax_ind_q(2) = sum(dmarg_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        dmarg_tax_ind_q(3) = sum(dmarg_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        dmarg_tax_ind_q(4) = sum(dmarg_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        dmarg_tax_ind_q(5) = sum(dmarg_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute average transfer rates by quintile (CBO concept: total transfers of quintile divided by total income of quintile)
        avg_transfers_q(1) = sum(transfers_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_transfers_q(2) = sum(transfers_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_transfers_q(3) = sum(transfers_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_transfers_q(4) = sum(transfers_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_transfers_q(5) = sum(transfers_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        ! compute average labor income tax rates by quintile (CBO concept: total labor tax of quintile divided by total income of quintile)
        avg_income_tax_q(1) = sum(income_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_income_tax_q(2) = sum(income_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_income_tax_q(3) = sum(income_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_income_tax_q(4) = sum(income_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_income_tax_q(5) = sum(income_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))


        ! compute average labor income tax rates by quintile (individual concept)
        avg_income_tax_ind_q(1) = sum((income_tax_sorted(1:p_20_ind)/ytot_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_income_tax_ind_q(2) = sum((income_tax_sorted(p_20_ind+1:p_40_ind)/ytot_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_income_tax_ind_q(3) = sum((income_tax_sorted(p_40_ind+1:p_60_ind)/ytot_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_income_tax_ind_q(4) = sum((income_tax_sorted(p_60_ind+1:p_80_ind)/ytot_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_income_tax_ind_q(5) = sum((income_tax_sorted(p_80_ind+1:Nstate)/ytot_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        ! compute average transfer rates by quintile (individual concept)
        avg_transfer_ind_q(1) = sum((transfers_sorted(1:p_20_ind)/ytot_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_transfer_ind_q(2) = sum((transfers_sorted(p_20_ind+1:p_40_ind)/ytot_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_transfer_ind_q(3) = sum((transfers_sorted(p_40_ind+1:p_60_ind)/ytot_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_transfer_ind_q(4) = sum((transfers_sorted(p_60_ind+1:p_80_ind)/ytot_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_transfer_ind_q(5) = sum((transfers_sorted(p_80_ind+1:Nstate)/ytot_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))


        ! compute average capital income tax rates by quintile (CBO concept: total capital tax of quintile divided by total income of quintile)
        avg_cap_income_tax_q(1) = sum(capital_income_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_cap_income_tax_q(2) = sum(capital_income_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_cap_income_tax_q(3) = sum(capital_income_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_cap_income_tax_q(4) = sum(capital_income_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_cap_income_tax_q(5) = sum(capital_income_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        ! compute average total tax rates (without transfers) by quintile (CBO concept: total tax of quintile divided by total income of quintile)
        avg_total_tax_q(1) = sum(total_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_total_tax_q(2) = sum(total_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_total_tax_q(3) = sum(total_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_total_tax_q(4) = sum(total_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_total_tax_q(5) = sum(total_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, 'AVERAGE TAX AND TRANSFER RATES (BY TOTAL INCOME)'
        print *, '** Average transfers rates by income quintile in percent, CBO concept =', avg_transfers_q*100d0
        print *, '** Average transfers rates by income quintile in percent; individual =', avg_transfer_ind_q*100d0
        print *, '** Average tax rates (without transfers) by income quintile in percent =', avg_total_tax_q*100d0
        print *, '** Average labor income tax rates by income quintile, CBO concept =', avg_income_tax_q*100d0
        print *, '** Average labor income tax rates by income quintile, individual =', avg_income_tax_ind_q*100d0
        print *, '** Average capital income tax rates by income quintile =', avg_cap_income_tax_q*100d0

        print *, '**'
        print *, '** Average rates,  CBO, total system = ', avg_tax_q*100D0
        print *, '** Average rates,  ind, total system = ', avg_tax_ind_q*100D0
        print *, '** Marginal rates, ind, total system = ', marg_tax_ind_q*100D0
        print *, '** DMarginal rates, ind, total system = ', dmarg_tax_ind_q*100D0
        print *, '**'



    
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Income Distribution
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        print *, 'INCOME AND WEALTH DISTRIBUTIONS'
        print *, ''
    
        ! Total income shares (CBO)
        p_90_ind = minloc(abs(mu_cumu-0.9d0),dim=1)
        p_99_ind = minloc(abs(mu_cumu-0.99d0),dim=1)

        inc_share(1) = sum(ytot_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(ytot_sorted*mu_sorted)
        inc_share(2) = sum(ytot_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(ytot_sorted*mu_sorted)
        inc_share(3) = sum(ytot_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(ytot_sorted*mu_sorted)
        inc_share(4) = sum(ytot_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(ytot_sorted*mu_sorted)
        inc_share(5) = sum(ytot_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(ytot_sorted*mu_sorted)
        inc_share(6) = sum(ytot_sorted(p_90_ind+1:Nstate)*mu_sorted(p_90_ind+1:Nstate))/sum(ytot_sorted*mu_sorted)
        inc_share(7)= sum(ytot_sorted(p_99_ind+1:Nstate)*mu_sorted(p_99_ind+1:Nstate))/sum(ytot_sorted*mu_sorted)
        
        print *, ' ** DIST TOT INC CBO **'
        print *, ' inc share = ', inc_share*100d0
        
        ! Labor income shares by total income groups
        yl_sorted = yl(ind_sort)
        yl_share(1) = sum(yl_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(yl_sorted*mu_sorted)
        yl_share(2) = sum(yl_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(yl_sorted*mu_sorted)
        yl_share(3) = sum(yl_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(yl_sorted*mu_sorted)
        yl_share(4) = sum(yl_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(yl_sorted*mu_sorted)
        yl_share(5) = sum(yl_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(yl_sorted*mu_sorted)
        
        print *, ''
        print *, 'Labor income shares for 5 quintiles of total income'
        print *, '** Labor income shares in percent = ', yl_share*100d0
    
        ! Capital income shares by total income groups
        yk_sorted = yk(ind_sort)
        yk_share(1) = sum(yk_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(yk_sorted*mu_sorted)
        yk_share(2) = sum(yk_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(yk_sorted*mu_sorted)
        yk_share(3) = sum(yk_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(yk_sorted*mu_sorted)
        yk_share(4) = sum(yk_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(yk_sorted*mu_sorted)
        yk_share(5) = sum(yk_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(yk_sorted*mu_sorted)
        
        print *, ''
        print *, 'Capital income shares for 5 quintiles of total income'
        print *, '** Capital income shares in percent = ', yk_share*100d0

        ! Labor income distributon
        ! Sort labor income
        yl_sorted = yl
        ind_sort = 0
        call hpsort_eps_epw(Nstate, yl_sorted, ind_sort, 1D-16)
        ! Sort measure correspondingly
        mu_sorted = muu(ind_sort)
        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do
        ! Indices of percentiles of the labor income distribution
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        p_99_ind = minloc(abs(mu_cumu-0.99d0),dim=1)
        ! More percentiles
        p_1_ind = minloc(abs(mu_cumu-0.01d0),dim=1)
        p_5_ind = minloc(abs(mu_cumu-0.05d0),dim=1)
        p_10_ind = minloc(abs(mu_cumu-0.10d0),dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.50d0),dim=1)

        ! Income shares of quintiles
        inc_share(1) = sum(yl_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(yl_sorted*mu_sorted)
        inc_share(2) = sum(yl_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(yl_sorted*mu_sorted)
        inc_share(3) = sum(yl_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(yl_sorted*mu_sorted)
        inc_share(4) = sum(yl_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(yl_sorted*mu_sorted)
        inc_share(5) = sum(yl_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(yl_sorted*mu_sorted)
        inc_share(6) = sum(yl_sorted(p_90_ind+1:Nstate)*mu_sorted(p_90_ind+1:Nstate))/sum(yl_sorted*mu_sorted)
        inc_share(7)= sum(yl_sorted(p_99_ind+1:Nstate)*mu_sorted(p_99_ind+1:Nstate))/sum(yl_sorted*mu_sorted)

        print *, ''
        print *, 'Labor income distribution for 5 quintiles, top 10%, and top 1%'
        print *, '** Labor income shares in percent = ', inc_share*100d0
        
        print *, ''
        print *, '** 1st percentile labor income rel. to median = ', yl_sorted(p_1_ind)/yl_sorted(p_50_ind)
        print *, '** 5th percentile labor income rel. to median = ', yl_sorted(p_5_ind)/yl_sorted(p_50_ind)
        print *, '** 10th percentile labor income rel. to median = ', yl_sorted(p_10_ind)/yl_sorted(p_50_ind)
        print *, '** 20th percentile labor income rel. to median = ', yl_sorted(p_20_ind)/yl_sorted(p_50_ind)
        
        ! compute average labor income tax rates by quintile: here labor income quintiles
        income_tax_sorted = income_tax(ind_sort)
        avg_income_tax_q(1) = sum(income_tax_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(yl_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))
        avg_income_tax_q(2) = sum(income_tax_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(yl_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))
        avg_income_tax_q(3) = sum(income_tax_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(yl_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))
        avg_income_tax_q(4) = sum(income_tax_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(yl_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))
        avg_income_tax_q(5) = sum(income_tax_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(yl_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, 'Labor income tax rates for labor income quintiles'
        print *, '** Tax rates in % = ', avg_income_tax_q*100d0


        av_tax_sorted = av_tax(ind_sort)                            ! average tax-and-transfer rate sorted by total income
        marg_tax_sorted = marg_tax(ind_sort)                        ! marginal tax-and-transfer rate sorted by total income
        tax_paid_sorted = tax_paid(ind_sort)
        
        print *, ''
        print *, '** Marginal t&T'
        print *, '** Marginal top 10 = ', sum(marg_tax_sorted(p_90_ind+1:Nstate)*mu_sorted(p_90_ind+1:Nstate))/sum(mu_sorted(p_90_ind+1:Nstate))
        print *, '** Marginal top 1 = ', sum(marg_tax_sorted(p_99_ind+1:Nstate)*mu_sorted(p_99_ind+1:Nstate))/sum(mu_sorted(p_99_ind+1:Nstate))
        print *, '** Marginal overall = ', sum(marg_tax_sorted*mu_sorted)/sum(mu_sorted)
        print *, '** Marginal top 10 - Marginal overall = ', sum(marg_tax_sorted(p_90_ind+1:Nstate)*mu_sorted(p_90_ind+1:Nstate))/sum(mu_sorted(p_90_ind+1:Nstate)) - sum(marg_tax_sorted*mu_sorted)/sum(mu_sorted)
        

        print *, '** Average labor tax'
        print *, '** Average overall (CBO) = ', sum(income_tax_sorted*mu_sorted)/sum(yl_sorted*mu_sorted)
        print *, '** Average overall (ind) = ', sum((income_tax_sorted/yl_sorted)*mu_sorted)/sum(mu_sorted)

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Net Worth Distribution
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
        ! Sort by assets
        a_sorted = S1u
        ind_sort = 0
        call hpsort_eps_epw(Nstate, a_sorted, ind_sort, 1D-16)
        ! sort measure correspondingly
        mu_sorted = muu(ind_sort)
        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do
        ! Compute indices of percentiles of the net worth distribution
        p_20_ind = minloc(abs(mu_cumu-0.2d0),dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0),dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0),dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0),dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0),dim=1)
        p_99_ind = minloc(abs(mu_cumu-0.99d0),dim=1)
    
        ! Net worth shares of quintiles, top 10%, and top 1%
        a_share(1) = sum(a_sorted(1:p_20_ind)*mu_sorted(1:p_20_ind))/sum(a_sorted*mu_sorted)
        a_share(2) = sum(a_sorted(p_20_ind+1:p_40_ind)*mu_sorted(p_20_ind+1:p_40_ind))/sum(a_sorted*mu_sorted)
        a_share(3) = sum(a_sorted(p_40_ind+1:p_60_ind)*mu_sorted(p_40_ind+1:p_60_ind))/sum(a_sorted*mu_sorted)
        a_share(4) = sum(a_sorted(p_60_ind+1:p_80_ind)*mu_sorted(p_60_ind+1:p_80_ind))/sum(a_sorted*mu_sorted)
        a_share(5) = sum(a_sorted(p_80_ind+1:Nstate)*mu_sorted(p_80_ind+1:Nstate))/sum(a_sorted*mu_sorted)
        a_share(6) = sum(a_sorted(p_90_ind+1:Nstate)*mu_sorted(p_90_ind+1:Nstate))/sum(a_sorted*mu_sorted)
        a_share(7) = sum(a_sorted(p_99_ind+1:Nstate)*mu_sorted(p_99_ind+1:Nstate))/sum(a_sorted*mu_sorted)
    
        print *, ''
        print *, 'Net worth distribution for 5 quintiles, top 10%, and top 1%'
        print *, '** Net worth shares in percent = ', a_share*100d0
    
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Mean and median capital income
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
        ! Sort by capital income
        yk_sorted = yk
        ind_sort = 0
        call hpsort_eps_epw(Nstate, yk_sorted, ind_sort, 1D-16)
        ! sort measure accordingly
        mu_sorted = muu(ind_sort)
        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do
        ! find index for median
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        !print *, ''
        !print *, 'Mean and median capital income relative to calibrated median income'
        !print *, 'Median capital income (rescaled by median total income) = ', yk_sorted(p_50_ind)/ymean
        !print *, 'Mean capital income (rescaled by median total income) = ', sum(yk*muu)/ymean


    end subroutine income_distribution
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution of tax-and-transfer rates for plots
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine implied_rates(index_ss,ngrid)
    
        implicit none
        
        ! Local variable declarations
        integer, intent(in) :: index_ss                 ! =1 if initial, 2 if LR
        integer, intent(in) :: ngrid                    ! number of points on which rates are computed
        real(8) :: yl_min, yl_max                       ! bounds for labor income grid
        real(8), dimension(ngrid) :: yl, yl_scaled      ! labor income grid, labor income grid relative to ymean
        real(8), dimension(ngrid) :: avg_tax_transfer   ! average tax-and-transfer rate
        real(8), dimension(ngrid) :: avg_tax            ! average tax rate
        real(8), dimension(ngrid) :: avg_transfer       ! average transfer
        real(8), dimension(ngrid) :: marg_tax_transfer  ! marginal tax-and-transfer rate
        real(8), dimension(ngrid) :: marg_tax           ! marginal tax rate
        real(8), dimension(ngrid) :: marg_transfer      ! marginal transfer rate
        integer :: iter                                 ! counter
        real(8) :: yk_mean                              ! mean capital income
        real(8), dimension(ngrid) :: et                 ! auxiliary variable transfer
        real(8) :: chit                                 ! auxiliary variable transfer
        
        ! Generate labor income grid
        yl_min = 0.001d0                        ! minimum for labor income grid
        yl_max = 5d0                            ! maximum for labor income grid
        yl = LINSPACE_FET(yl_min,yl_max,ngrid)  ! labor income grid
        yl_scaled = yl*ymean                    ! rescale labor income grid
        
        open(101,  file = 'Output/yl_grid.txt',      status = 'unknown')
        write(101, *) yl
        close(101)
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Case 1: zero capital income
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        ! Total tax-and-transfer rates
        do iter = 1, ngrid
            avg_tax_transfer(iter) = TAX_sc(yl_scaled(iter),0d0)/yl_scaled(iter)
            marg_tax_transfer(iter) = dTAX_sc(yl_scaled(iter),0d0)
        enddo
        
        ! Labor tax rates
        avg_tax = (exp(log(lambda)*((yl_scaled/ymean)**(-gamma)))*yl_scaled)/yl_scaled
        marg_tax = exp(log(lambda)*((yl_scaled/ymean)**(-gamma)))*(1d0-gamma*log(lambda)*((yl_scaled/ymean)**(-gamma)))
        
        ! Transfer rates
        et   = exp(-rt*((yl_scaled)/ymean-omegat))                          
        chit = exp(rt*omegat)/(1.0D0 + exp(rt*omegat))                      
        avg_transfer = -((mt*ymean)*(et/(1.0D0+et))*(1.0D0/chit))/yl_scaled  
        marg_transfer = (rt/ymean)*(mt*ymean)*(et/((1.0D0+et)**2d0))*(1D0/chit)
        
        if (index_ss == 1) then
            
            open(102,  file = 'Output/avg_tax_transfer_zerok_init.txt',      status = 'unknown')
            write(102, *) avg_tax_transfer
            close(102)
            
            open(103,  file = 'Output/marg_tax_transfer_zerok_init.txt',      status = 'unknown')
            write(103, *) marg_tax_transfer
            close(103)
            
            open(104,  file = 'Output/avg_tax_zerok_init.txt',      status = 'unknown')
            write(104, *) avg_tax
            close(104)
            
            open(105,  file = 'Output/marg_tax_zerok_init.txt',      status = 'unknown')
            write(105, *) marg_tax
            close(105)
            
            open(106,  file = 'Output/avg_transfer_zerok_init.txt',      status = 'unknown')
            write(106, *) avg_transfer
            close(106)
            
            open(107,  file = 'Output/marg_transfer_zerok_init.txt',      status = 'unknown')
            write(107, *) marg_transfer
            close(107)
            
        else
            
            open(102,  file = 'Output/avg_tax_transfer_zerok_lr.txt',      status = 'unknown')
            write(102, *) avg_tax_transfer
            close(102)
            
            open(103,  file = 'Output/marg_tax_transfer_zerok_lr.txt',      status = 'unknown')
            write(103, *) marg_tax_transfer
            close(103)
            
            open(104,  file = 'Output/avg_tax_zerok_lr.txt',      status = 'unknown')
            write(104, *) avg_tax
            close(104)
            
            open(105,  file = 'Output/marg_tax_zerok_lr.txt',      status = 'unknown')
            write(105, *) marg_tax
            close(105)
            
            open(106,  file = 'Output/avg_transfer_zerok_lr.txt',      status = 'unknown')
            write(106, *) avg_transfer
            close(106)
            
            open(107,  file = 'Output/marg_transfer_zerok_lr.txt',      status = 'unknown')
            write(107, *) marg_transfer
            close(107)
            
        endif
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Case 2: mean capital income
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        ! Total tax-and-transfer rates
        yk_mean = sum(mu*r*S(:,1))
        do iter = 1, ngrid
            avg_tax_transfer(iter) = TAX_sc(yl_scaled(iter),yk_mean)/yl_scaled(iter)
            marg_tax_transfer(iter) = dTAX_sc(yl_scaled(iter),yk_mean)
        enddo
        
        ! Labor tax rates
        avg_tax = (exp(log(lambda)*((yl_scaled/ymean)**(-gamma)))*yl_scaled)/yl_scaled
        marg_tax = exp(log(lambda)*((yl_scaled/ymean)**(-gamma)))*(1d0-gamma*log(lambda)*((yl_scaled/ymean)**(-gamma)))
        
        ! Transfer rates
        et   = exp(-rt*((yl_scaled+yk_mean)/ymean-omegat))                          
        chit = exp(rt*omegat)/(1.0D0 + exp(rt*omegat))                      
        avg_transfer = -((mt*ymean)*(et/(1.0D0+et))*(1.0D0/chit))/yl_scaled  
        marg_transfer = (rt/ymean)*(mt*ymean)*(et/((1.0D0+et)**2d0))*(1D0/chit)
        
        if (index_ss == 1) then
            
            open(102,  file = 'Output/avg_tax_transfer_meank_init.txt',      status = 'unknown')
            write(102, *) avg_tax_transfer
            close(102)
            
            open(103,  file = 'Output/marg_tax_transfer_meank_init.txt',      status = 'unknown')
            write(103, *) marg_tax_transfer
            close(103)
            
            open(104,  file = 'Output/avg_tax_meank_init.txt',      status = 'unknown')
            write(104, *) avg_tax
            close(104)
            
            open(105,  file = 'Output/marg_tax_meank_init.txt',      status = 'unknown')
            write(105, *) marg_tax
            close(105)
            
            open(106,  file = 'Output/avg_transfer_meank_init.txt',      status = 'unknown')
            write(106, *) avg_transfer
            close(106)
            
            open(107,  file = 'Output/marg_transfer_meank_init.txt',      status = 'unknown')
            write(107, *) marg_transfer
            close(107)
            
        else
            
            open(102,  file = 'Output/avg_tax_transfer_meank_lr.txt',      status = 'unknown')
            write(102, *) avg_tax_transfer
            close(102)
            
            open(103,  file = 'Output/marg_tax_transfer_meank_lr.txt',      status = 'unknown')
            write(103, *) marg_tax_transfer
            close(103)
            
            open(104,  file = 'Output/avg_tax_meank_lr.txt',      status = 'unknown')
            write(104, *) avg_tax
            close(104)
            
            open(105,  file = 'Output/marg_tax_meank_lr.txt',      status = 'unknown')
            write(105, *) marg_tax
            close(105)
            
            open(106,  file = 'Output/avg_transfer_meank_lr.txt',      status = 'unknown')
            write(106, *) avg_transfer
            close(106)
            
            open(107,  file = 'Output/marg_transfer_meank_lr.txt',      status = 'unknown')
            write(107, *) marg_transfer
            close(107)
            
        endif
        
        
    end subroutine implied_rates
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution marginal tax rates
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    subroutine tax_distribution(index_ss)
    
        implicit none
        
        ! Local variable declarations
        integer, intent(in) :: index_ss                         ! =1 if initial, 2 if LR
        real(8), dimension(NS) :: yl, yk, ytot                  ! labor, capital, total income
        real(8), dimension(NS) :: marg_rate_dist, avg_rate_dist ! distribution of total marginal and average t&T rates
        integer :: is                                           ! counter   
        real(8), dimension(NS) :: marg_tax_dist, avg_tax_dist   ! distribution of total marginal and average tax rates
        real(8), dimension(NS) :: marg_transfer_dist, avg_transfer_dist   ! distribution of total marginal and average transfer rates
        real(8), dimension(NS) :: et                            ! auxiliary variable transfer
        real(8) :: chit                                         ! auxiliary variable transfer
        
        ! Compute tax-transfer rates across distribution
        yl = wge*h_pol*S(:,2)       ! labor income
        yk = r*S(:,1)               ! capital income
        ytot = yl + yk              ! total income
        do is = 1, NS
            avg_rate_dist(is) = TAX_sc(yl(is),yk(is))/ytot(is)
            marg_rate_dist(is) = dTAX_sc(yl(is),yk(is))
        enddo
        
        ! Labor tax rates
        avg_tax_dist = (exp(log(lambda)*((yl/ymean)**(-gamma)))*yl)/ytot
        marg_tax_dist = exp(log(lambda)*((yl/ymean)**(-gamma)))*(1d0-gamma*log(lambda)*((yl/ymean)**(-gamma)))
        
        ! Transfer rates
        et   = exp(-rt*((yl+yk)/ymean-omegat))                          
        chit = exp(rt*omegat)/(1.0D0 + exp(rt*omegat))                      
        avg_transfer_dist = -((mt*ymean)*(et/(1.0D0+et))*(1.0D0/chit))/ytot  
        marg_transfer_dist = (rt/ymean)*(mt*ymean)*(et/((1.0D0+et)**2d0))*(1D0/chit)
        
        ! Save distributions
        if (index_ss == 1) then
            
            open(101,  file = 'Output/avg_rate_dist_init.txt',      status = 'unknown')
            write(101, *) avg_rate_dist
            close(101)
            
            open(102,  file = 'Output/marg_rate_dist_init.txt',      status = 'unknown')
            write(102, *) marg_rate_dist
            close(102)
            
            open(103,  file = 'Output/avg_tax_dist_init.txt',      status = 'unknown')
            write(103, *) avg_tax_dist
            close(103)
            
            open(104,  file = 'Output/marg_tax_dist_init.txt',      status = 'unknown')
            write(104, *) marg_tax_dist
            close(104)
            
            open(105,  file = 'Output/avg_transfer_dist_init.txt',      status = 'unknown')
            write(105, *) avg_transfer_dist
            close(105)
            
            open(106,  file = 'Output/marg_transfer_dist_init.txt',      status = 'unknown')
            write(106, *) marg_transfer_dist
            close(106)
            
        else
            
            open(101,  file = 'Output/avg_rate_dist_lr.txt',      status = 'unknown')
            write(101, *) avg_rate_dist
            close(101)
            
            open(102,  file = 'Output/marg_rate_dist_lr.txt',      status = 'unknown')
            write(102, *) marg_rate_dist
            close(102)
            
            open(103,  file = 'Output/avg_tax_dist_lr.txt',      status = 'unknown')
            write(103, *) avg_tax_dist
            close(103)
            
            open(104,  file = 'Output/marg_tax_dist_lr.txt',      status = 'unknown')
            write(104, *) marg_tax_dist
            close(104)
            
            open(105,  file = 'Output/avg_transfer_dist_lr.txt',      status = 'unknown')
            write(105, *) avg_transfer_dist
            close(105)
            
            open(106,  file = 'Output/marg_transfer_dist_lr.txt',      status = 'unknown')
            write(106, *) marg_transfer_dist
            close(106)
        
        endif
        
    end subroutine tax_distribution

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution labor earnings growth rates
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine income_growth

        implicit none

        ! Local variable declarations
        real(8) :: inc(NS*Nx)           ! today's income for every state (#NS) and possible productivity tomorrow (#Nx) [they all choose the same a]
        real(8) :: inc_p(NS*Nx)         ! tomorrow's income for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: inc_d(NS*Nx)         ! income difference for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: inc_d_sorted(NS*Nx)  ! sorted income difference for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: mass(NS*Nx)          ! mass for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: mass_sorted(NS*Nx)   ! sorted mass
        real(8) :: mass_cum(NS*Nx)      ! cumulative mass
        integer :: ind_sort(NS*Nx)      ! sorting indicator
        integer :: iter, iter_s, iter_xp, is_np1, is_np2    ! counter
        real(8) :: inc_aux, mu_aux, ap_aux, ap_aux_vec(1), da, h_aux    ! auxiliary real variables
        integer :: ap_ind(2), ap_ind_aux(1)                 ! auxiliary integer variables
        real(8) :: prob_x               ! auxiliary real variable
        integer :: p_8_ind, p_10_ind, p_12_ind, p_20_ind, p_30_ind, p_40_ind, p_50_ind, p_60_ind, p_70_ind, p_80_ind, p_90_ind, p_25_ind, p_75_ind, p_2_5_ind, p_97_5_ind  ! indices for percentiles
        real(8) :: m1, m2, m3, m4       ! moments
        real(8) :: inc_d_mean           ! mean of earnings growth distribution
        real(8) :: inc_d_var            ! variance
        real(8) :: inc_d_std            ! standard deviation
        real(8) :: inc_d_skew           ! skewness
        real(8) :: inc_d_kurt           ! kurtosis
        real(8) :: inc_d_9010           ! P90-P10
        real(8) :: inc_d_kelley         ! Kelley skewness
        real(8) :: inc_d_crow           ! Crow-Siddiqui kurtosis
        real(8) :: log_inc(NS*Nx)       ! log of today's income for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: log_inc_var          ! variance of log income
        real(8) :: log_inc_sd           ! standard deviation of log income
        real(8) :: m1_log_inc, m2_log_inc   ! moments 

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Compute distribution of earnings growth rates
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        ! initialize counter
        iter = 0                
        
        ! loop over all states today
        do iter_s = 1, NS       
            
            ! Current income
            inc_aux = wge*h_pol(iter_s)*x_vec(xind(iter_s))
            ! Mass at this point
            mu_aux = mu(iter_s)
            ! Asset choice
            ap_aux = min(a_pol(iter_s),amax)
            ap_aux_vec = ap_aux
            ap_ind_aux = vsearchprm_FET(ap_aux_vec,a_vec,aprm)
            ap_ind(1) = ap_ind_aux(1)
            ap_ind(2) = ap_ind(1)+1
            da = (ap_aux - a_vec(ap_ind(1)))/(a_vec(ap_ind(1)+1) - a_vec(ap_ind(1)))

            ! Loop over all possible future productivities
            do iter_xp = 1, Nx
                            
                ! Add to iteration counter
                iter = iter+1

                ! Today's income
                inc(iter) = inc_aux
                
                ! Indices tomorrow given asset choice
                is_np1 = ap_ind(1) + Na*(iter_xp - 1)
                is_np2 = ap_ind(2) + Na*(iter_xp - 1)
                
                ! Probability of this x-transition
                prob_x = Px(xind(iter_s),iter_xp)
                
                ! Interpolate labor policy
                h_aux = (1d0-da)*h_pol(is_np1)+da*h_pol(is_np2)

                ! Future Income
                inc_p(iter) = wge*h_aux*x_vec(xind(is_np1))
                
                ! Mass in this state that makes this x-transition
                mass(iter) = mu_aux*prob_x

                ! Income growth for this state with this x-transition
                inc_d(iter) = log(inc_p(iter)) - log(inc(iter))
               
            ! End loop over potential future productivities
            enddo
        
        ! End loop over states
        enddo
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Compute moments of earnings growth distribution
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        ! Sort earnings 
        inc_d_sorted = inc_d
        ind_sort = 0
        call hpsort_eps_epw(NS*Nx,inc_d_sorted,ind_sort,1D-16)

        ! Sort mass correspondingly
        mass_sorted = mass(ind_sort)

        ! Compute cumulative mass
        mass_cum(1) = mass_sorted(1)
        do iter = 2, NS*Nx
            mass_cum(iter) = mass_cum(iter-1) + mass_sorted(iter)
        enddo

        ! Compute percentiles of labor earnings growth distribution
        p_8_ind = minloc(abs(mass_cum-0.08d0), dim=1)
        p_10_ind = minloc(abs(mass_cum-0.1d0), dim=1)
        p_12_ind = minloc(abs(mass_cum-0.12d0), dim=1)
        p_20_ind = minloc(abs(mass_cum-0.2d0), dim=1)
        p_30_ind = minloc(abs(mass_cum-0.3d0), dim=1)
        p_40_ind = minloc(abs(mass_cum-0.4d0), dim=1)
        p_50_ind = minloc(abs(mass_cum-0.5d0), dim=1)
        p_60_ind = minloc(abs(mass_cum-0.6d0), dim=1)
        p_70_ind = minloc(abs(mass_cum-0.7d0), dim=1)
        p_80_ind = minloc(abs(mass_cum-0.8d0), dim=1)
        p_90_ind = minloc(abs(mass_cum-0.9d0), dim=1)
        p_75_ind = minloc(abs(mass_cum-0.75d0), dim=1)
        p_25_ind = minloc(abs(mass_cum-0.25d0), dim=1)
        p_97_5_ind = minloc(abs(mass_cum-0.975d0), dim=1)
        p_2_5_ind = minloc(abs(mass_cum-0.025d0), dim=1)

        ! Compute moments of labor earnings growth distribution
        m1 = sum(inc_d_sorted*mass_sorted)
        m2 = sum(mass_sorted*(inc_d_sorted-m1)**2d0)
        m3 = sum(mass_sorted*(inc_d_sorted-m1)**3d0)
        m4 = sum(mass_sorted*(inc_d_sorted-m1)**4d0)
        inc_d_mean = m1
        inc_d_var = m2
        inc_d_std = sqrt(m2)
        inc_d_skew = m3*m2**(-3d0/2d0)
        inc_d_kurt = m4*m2**(-2d0)
        inc_d_9010 = inc_d_sorted(p_90_ind)-inc_d_sorted(p_10_ind)
        inc_d_kelley = ((inc_d_sorted(p_90_ind)-inc_d_sorted(p_50_ind))-(inc_d_sorted(p_50_ind)-inc_d_sorted(p_10_ind)))/(inc_d_sorted(p_90_ind)-inc_d_sorted(p_10_ind))
        inc_d_crow = (inc_d_sorted(p_97_5_ind)-inc_d_sorted(p_2_5_ind))/(inc_d_sorted(p_75_ind)-inc_d_sorted(p_25_ind))

        print *, ''
        print *, "Earnings Growth Distribution"
        print *, "Mean = ", inc_d_mean
        print *, "Variance = ", inc_d_var
        print *, "P90-P10 = ", inc_d_9010
        print *, "Std. Dev. = ", inc_d_std
        print *, "Skewness = ", inc_d_skew
        print *, "Kelley = ", inc_d_kelley
        print *, "Kurtosis = ", inc_d_kurt
        print *, "Crow Siddiqui = ", inc_d_crow
        print *, "10th percentile = ", inc_d_sorted(p_10_ind)
        !print *, "8-12th percentile = ", sum(mass_sorted(p_8_ind:p_12_ind)*inc_d_sorted(p_8_ind:p_12_ind))/sum(mass_sorted(p_8_ind:p_12_ind))
        print *, "20th percentile = ", inc_d_sorted(p_20_ind)
        print *, "30th percentile = ", inc_d_sorted(p_30_ind)
        print *, "40th percentile = ", inc_d_sorted(p_40_ind)
        print *, "50th percentile = ", inc_d_sorted(p_50_ind)
        print *, "60th percentile = ", inc_d_sorted(p_60_ind)
        print *, "70th percentile = ", inc_d_sorted(p_70_ind)
        print *, "80th percentile = ", inc_d_sorted(p_80_ind)
        print *, "90th percentile = ", inc_d_sorted(p_90_ind)

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Variance of log income
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        log_inc = log(inc)
        m1_log_inc = sum(log_inc*mass)
        m2_log_inc = sum(mass*(log_inc-m1_log_inc)**2d0)
        log_inc_var = m2_log_inc
        log_inc_sd = sqrt(m2_log_inc)
        print *, ''
        print *, "Variance of log income = ", log_inc_var

    end subroutine income_growth


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Distribution hours
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine hours_distribution

        implicit none

        ! Local variable declarations
        real(8), dimension(NS) :: yl, ytot, yh                                                  ! labor income, total income, hourly wage
        real(8), dimension(NS) :: ytot_sorted, yl_sorted, yh_sorted     ! total income, labor income, tax paid, transfers (all sorted)
        real(8), dimension(NS) :: mu_sorted, mu_cumu, a_sorted     ! measure sorted, cum. measure, assets sorted, cap.+lab. tax, also sorted
        real(8), dimension(NS) :: h_sorted, h_normalized
        integer :: ind_sort(NS)                                                                 ! index for sorting
        integer :: Nstate                                                                       ! auxiliary variable equal to NS
        real(8), dimension(NS) :: muu, S1u, S2u                                                 ! auxiliary variables equal to measure, asset states, productivity states
        real(8) :: avg_h_dec(10), avg_h_q(5)
        integer :: p_20_ind, p_40_ind, p_60_ind, p_80_ind, p_90_ind                              ! indices for percentiles
        integer :: p_10_ind, p_30_ind, p_50_ind, p_70_ind                                      ! indices for percentiles
        integer :: iter, is                                                                     ! counter

        real(8) :: log_h(NS)       ! log of today's income for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: log_h_var          ! variance of log income
        real(8) :: log_h_sd           ! standard deviation of log income
        real(8) :: m1_log_h, m2_log_h   ! moments

        real(8) :: log_c(NS)            ! log of today's income for every state (#NS) and possible productivity tomorrow (#Nx)
        real(8) :: log_c_var            ! variance of log income
        real(8) :: log_c_sd             ! standard deviation of log income
        real(8) :: m1_log_c, m2_log_c   ! moments

        real(8) :: aux_x, aux_y, corr_hw, mean_hours_lev_q(5)

        ! auxiliary variables
        Nstate = NS     ! number of individual states
        muu = mu        ! measure over individual states
        S1u = S(:,1)    ! assets for individual states
        S2u = S(:,2)    ! productivities for individual states

        ! Income variables
        yh = wge*S2u        ! labor income
        
        
        ! 
        aux_x = sum((h_pol-sum(mu*h_pol))*(S2u-sum(mu*S2u))*mu)
        aux_y = sqrt(sum(((h_pol-sum(mu*h_pol))**2d0)*mu))*sqrt(sum(((S2u-sum(mu*S2u))**2d0)*mu))
        corr_hw = aux_x/aux_y
        print *, "corr_hw = ", corr_hw
        

        ! Normalize hours by median hours
        h_sorted = h_pol
        ind_sort = 0
        call hpsort_eps_epw(Nstate, h_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_30_ind = minloc(abs(mu_cumu-0.3d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_70_ind = minloc(abs(mu_cumu-0.7d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        h_normalized = log(h_pol) - sum(mu*log(h_pol))!/Hagg!h_sorted(p_50_ind)


        ! mean hours by quintile
        mean_hours_lev_q(1) = sum((h_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        mean_hours_lev_q(2) = sum((h_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        mean_hours_lev_q(3) = sum((h_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        mean_hours_lev_q(4) = sum((h_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        mean_hours_lev_q(5) = sum((h_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))
        print *, "Mean hours in levels by quintile = ", mean_hours_lev_q

        ! sort hourly income array
        yh_sorted = yh
        ind_sort = 0
        call hpsort_eps_epw(Nstate, yh_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income
        h_sorted = h_normalized(ind_sort)

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_30_ind = minloc(abs(mu_cumu-0.3d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_70_ind = minloc(abs(mu_cumu-0.7d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        ! compute average hours  by decile
        avg_h_dec(1) = sum((h_sorted(1:p_10_ind))*mu_sorted(1:p_10_ind))/sum(mu_sorted(1:p_10_ind))
        avg_h_dec(2) = sum((h_sorted(p_10_ind+1:p_20_ind))*mu_sorted(p_10_ind+1:p_20_ind))/sum(mu_sorted(p_10_ind+1:p_20_ind))
        avg_h_dec(3) = sum((h_sorted(p_20_ind+1:p_30_ind))*mu_sorted(p_20_ind+1:p_30_ind))/sum(mu_sorted(p_20_ind+1:p_30_ind))
        avg_h_dec(4) = sum((h_sorted(p_30_ind+1:p_40_ind))*mu_sorted(p_30_ind+1:p_40_ind))/sum(mu_sorted(p_30_ind+1:p_40_ind))
        avg_h_dec(5) = sum((h_sorted(p_40_ind+1:p_50_ind))*mu_sorted(p_40_ind+1:p_50_ind))/sum(mu_sorted(p_40_ind+1:p_50_ind))
        avg_h_dec(6) = sum((h_sorted(p_50_ind+1:p_60_ind))*mu_sorted(p_50_ind+1:p_60_ind))/sum(mu_sorted(p_50_ind+1:p_60_ind))
        avg_h_dec(7) = sum((h_sorted(p_60_ind+1:p_70_ind))*mu_sorted(p_60_ind+1:p_70_ind))/sum(mu_sorted(p_60_ind+1:p_70_ind))
        avg_h_dec(8) = sum((h_sorted(p_70_ind+1:p_80_ind))*mu_sorted(p_70_ind+1:p_80_ind))/sum(mu_sorted(p_70_ind+1:p_80_ind))
        avg_h_dec(9) = sum((h_sorted(p_80_ind+1:p_90_ind))*mu_sorted(p_80_ind+1:p_90_ind))/sum(mu_sorted(p_80_ind+1:p_90_ind))
        avg_h_dec(10) = sum((h_sorted(p_90_ind+1:Nstate))*mu_sorted(p_90_ind+1:Nstate))/sum(mu_sorted(p_90_ind+1:Nstate))

        ! By quintile
        ! compute average hours  by decile
        avg_h_q(1) = sum((h_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_h_q(2) = sum((h_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_h_q(3) = sum((h_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_h_q(4) = sum((h_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_h_q(5) = sum((h_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, 'log(Hours/mean hours) by wage'
        !print *, '** by decile ', avg_h_dec
        print *, '** by quintile ', avg_h_q


        ! sort by hours
        h_sorted = h_normalized
        ind_sort = 0
        call hpsort_eps_epw(Nstate, h_sorted, ind_sort, 1D-16)
        ! sort measure and taxes paid arrays correspondingly
        mu_sorted = muu(ind_sort)                                   ! measure sorted by total income

        ! compute cumulative measure
        mu_cumu = 0d0
        mu_cumu(1) = mu_sorted(1)
        do iter = 2, Nstate
            mu_cumu(iter) = mu_cumu(iter-1) + mu_sorted(iter)
        end do

        p_10_ind = minloc(abs(mu_cumu-0.1d0), dim=1)
        p_20_ind = minloc(abs(mu_cumu-0.2d0), dim=1)
        p_30_ind = minloc(abs(mu_cumu-0.3d0), dim=1)
        p_40_ind = minloc(abs(mu_cumu-0.4d0), dim=1)
        p_50_ind = minloc(abs(mu_cumu-0.5d0), dim=1)
        p_60_ind = minloc(abs(mu_cumu-0.6d0), dim=1)
        p_70_ind = minloc(abs(mu_cumu-0.7d0), dim=1)
        p_80_ind = minloc(abs(mu_cumu-0.8d0), dim=1)
        p_90_ind = minloc(abs(mu_cumu-0.9d0), dim=1)
        !print *, 'Top quintile cutoff = ', ytot_sorted(p_80_ind)/ymean

        ! compute average hours  by decile
        avg_h_dec(1) = sum((h_sorted(1:p_10_ind))*mu_sorted(1:p_10_ind))/sum(mu_sorted(1:p_10_ind))
        avg_h_dec(2) = sum((h_sorted(p_10_ind+1:p_20_ind))*mu_sorted(p_10_ind+1:p_20_ind))/sum(mu_sorted(p_10_ind+1:p_20_ind))
        avg_h_dec(3) = sum((h_sorted(p_20_ind+1:p_30_ind))*mu_sorted(p_20_ind+1:p_30_ind))/sum(mu_sorted(p_20_ind+1:p_30_ind))
        avg_h_dec(4) = sum((h_sorted(p_30_ind+1:p_40_ind))*mu_sorted(p_30_ind+1:p_40_ind))/sum(mu_sorted(p_30_ind+1:p_40_ind))
        avg_h_dec(5) = sum((h_sorted(p_40_ind+1:p_50_ind))*mu_sorted(p_40_ind+1:p_50_ind))/sum(mu_sorted(p_40_ind+1:p_50_ind))
        avg_h_dec(6) = sum((h_sorted(p_50_ind+1:p_60_ind))*mu_sorted(p_50_ind+1:p_60_ind))/sum(mu_sorted(p_50_ind+1:p_60_ind))
        avg_h_dec(7) = sum((h_sorted(p_60_ind+1:p_70_ind))*mu_sorted(p_60_ind+1:p_70_ind))/sum(mu_sorted(p_60_ind+1:p_70_ind))
        avg_h_dec(8) = sum((h_sorted(p_70_ind+1:p_80_ind))*mu_sorted(p_70_ind+1:p_80_ind))/sum(mu_sorted(p_70_ind+1:p_80_ind))
        avg_h_dec(9) = sum((h_sorted(p_80_ind+1:p_90_ind))*mu_sorted(p_80_ind+1:p_90_ind))/sum(mu_sorted(p_80_ind+1:p_90_ind))
        avg_h_dec(10) = sum((h_sorted(p_90_ind+1:Nstate))*mu_sorted(p_90_ind+1:Nstate))/sum(mu_sorted(p_90_ind+1:Nstate))

        ! By quintile
        ! compute average hours  by decile
        avg_h_q(1) = sum((h_sorted(1:p_20_ind))*mu_sorted(1:p_20_ind))/sum(mu_sorted(1:p_20_ind))
        avg_h_q(2) = sum((h_sorted(p_20_ind+1:p_40_ind))*mu_sorted(p_20_ind+1:p_40_ind))/sum(mu_sorted(p_20_ind+1:p_40_ind))
        avg_h_q(3) = sum((h_sorted(p_40_ind+1:p_60_ind))*mu_sorted(p_40_ind+1:p_60_ind))/sum(mu_sorted(p_40_ind+1:p_60_ind))
        avg_h_q(4) = sum((h_sorted(p_60_ind+1:p_80_ind))*mu_sorted(p_60_ind+1:p_80_ind))/sum(mu_sorted(p_60_ind+1:p_80_ind))
        avg_h_q(5) = sum((h_sorted(p_80_ind+1:Nstate))*mu_sorted(p_80_ind+1:Nstate))/sum(mu_sorted(p_80_ind+1:Nstate))

        print *, ''
        print *, 'log(Hours/mean hours) by hours'
        !print *, '** by decile ', avg_h_dec
        print *, '** by quintile ', avg_h_q



        log_h = log(h_pol)
        m1_log_h = sum(log_h*muu)
        m2_log_h = sum(muu*(log_h-m1_log_h)**2d0)
        log_h_var = m2_log_h
        log_h_sd = sqrt(m2_log_h)
        print *, ''
        print *, "Variance of log h = ", log_h_var


        log_c = log(c_pol)
        m1_log_c = sum(log_c*muu)
        m2_log_c = sum(muu*(log_c-m1_log_c)**2d0)
        log_c_var = m2_log_c
        log_c_sd = sqrt(m2_log_c)
        print *, ''
        print *, "Variance of log c = ", log_c_var


    end subroutine hours_distribution

end module
